Distribution of the number of significant effect sizes

James E. Pustejovsky

A while back, I posted the outline of a problem about the number of significant effect size estimates in a study that reports multiple outcomes. This problem interests me because it connects to the issue of selective reporting of study results, which creates problems for meta-analysis. Here, I’ll re-state the problem in slightly more general terms and then make some notes about what’s going on.

Consider a study that assesses some effect size across $m$ different outcomes. (We’ll be thinking about one study at a time here, so no need to index the study as we would in a meta-analysis problem.) Let $T_i$ denote the effect size estimate for outcome $i$, let $V_i$ denote the sampling variance of the effect size estimate for outcome $i$, and let $\theta_i$ denote the true effect size parameter for corresponding to outcome $i$. Assume that the study outcomes $\left[T_i\right]_{i=1}^m$ follow a correlated-and-hierarchical effects model, in which \[T_i = \mu + u + v_i + e_i,\] where the study-level error $u \sim N\left(0,\tau^2\right)$, the effect-specific error $v_i \stackrel{iid}{\sim} N\left(0, \omega^2\right)$, and the vector of sampling errors $\left[e_i\right]_{i=1}^m$ is multivariate normal with mean $\mathbf{0}$, known variances $\text{Var}(e_i) = \sigma^2$, and compound symmetric correlation structure $\text{cor}(e_h, e_i) = \rho$.

Define $A_i$ as an indicator that is equal to one if $T_i$ is statistically significant at level $\alpha$ based on a one-sided test, and otherwise equal to zero. (Equivalently, let $A_i$ be equal to one if the effect is statistically significant at level $2 \alpha$ and in the theoretically expected direction.) Formally, \[A_i = I\left(\frac{T_i}{\sigma} > q_\alpha \right)\] where $q_\alpha = \Phi^{-1}(1 - \alpha)$ is the critical value from a standard normal distribution (e.g., $q_{.05} = 1.645$, $q_{.025} = 1.96$). Let $N_A = \sum_{i=1}^m A_i$ denote the total number of statistically significant effect sizes in the study. The question is: what is the distribution of $N_A$.

Compound symmetry to the rescue

As I noted in the previous post, this set-up means that the effect size estimates have a compound symmetric distribution. We can make this a bit more explicit by writing the sampling errors in terms of the sum of a component that’s common acrosss outcomes and a component that’s specific to each outcome. Thus, let $e_i = f + g_i$, where $f \sim N\left(0, \rho \sigma^2 \right)$ and $g_i \stackrel{iid}{\sim} N \left(0, (1 - \rho) \sigma^2\right)$. Let me also define $\zeta = \mu + u + f$ as the conditional mean of the effects. It then follows that the effect size estimates are conditionally independent, given the common components: \[ \left(T_i | \zeta \right) \stackrel{iid}{\sim} N\left(\zeta, \omega^2 + (1 - \rho) \sigma^2\right) \] Furthermore, the conditional probability of a significant effect is \[ \text{Pr}(A_i = 1 | \zeta) = \Phi\left(\frac{\zeta - q_{\alpha} \sigma}{\sqrt{\omega^2 + (1 - \rho)\sigma^2}}\right) \] and $A_1,...,A_m$ are mutually independent, conditional on $\zeta$. Therefore, the conditional distribution of $N_A$ is binomial, \[ \left(N_A | \zeta\right) \sim Bin(m, \pi) \] where \[ \pi = \Phi\left(\frac{\zeta - q_{\alpha} \sigma}{\sqrt{\omega^2 + (1 - \rho)\sigma^2}}\right). \] What about the unconditional distribution?

To get rid of the $\zeta$, we need to integrate over its distribution, which leads to \[ \text{Pr}(N_A = a) = \text{E}\left[\text{Pr}\left(N_A | \zeta\right)\right] = \int f_{N_A}\left(a | \zeta, \omega, \sigma, \rho, m\right) \times f_\zeta(\zeta | \mu, \tau, \sigma, \rho) \ d \zeta, \] where $f_{N_A}\left(a | \zeta, \omega, \sigma, \rho \right)$ is a binomial density with size $m$ and probability $\pi = \pi(\zeta, \omega, \sigma, \rho)$ and $f_\zeta(\zeta | \mu, \tau, \sigma, \rho)$ is a normal density with mean $\mu$ and variance $\tau^2 + \rho \sigma^2$.

This distribution is what you might call a binomial-normal convolution or a random-intercept probit model (where the random intercept is $\zeta$). As far as I know, the distribution cannot be evaluated analytically but instead must be calculated using some sort of numerical integration routine.

Just the moments, please

If all we care about is the expectation of $N_A$, we don’t need to bother with all the conditioning business and can just look at the marginal distribution of the effect size estimates taken individually. Marginally, $T_i$ is normally distributed with mean $\mu$ and variance $\tau^2 + \omega^2 + \sigma^2$, so $\text{Pr}(A_i = 1) = \psi$, where \[ \psi = \Phi\left(\frac{\mu - q_{\alpha} \sigma}{\sqrt{\tau^2 + \omega^2 + \sigma^2}}\right). \] By the linearity of expectations, \[ \text{E}(N_A) = \sum_{i=1}^m \text{E}(A_i) = m \psi. \]

We can also get an approximation for the variance of $N_A$ by working with its conditional distribution above. By the rule of variance decomposition, \[ \begin{aligned} \text{Var}(N_A) &= \text{E}\left[\text{Var}\left(N_A | \zeta\right)\right] + \text{Var}\left[\text{E}\left(N_A | \zeta\right)\right] \\ &= m \times \text{E}\left[\pi (1 - \pi)\right] + m^2 \times \text{Var}\left[\pi\right]\\ &= m \times \text{E}\left[\pi\right] \left(1 - \text{E}\left[\pi\right]\right) + m (m - 1) \times \text{Var}\left[\pi\right], \end{aligned} \] where $\pi$ is, as defined above, a function of $\zeta$ and thus a random variable. Now, $\text{E}(\pi) = \psi$ and we can get something close to $\text{Var}(\pi)$ using a first-order approximation: \[ \text{Var}\left(\pi\right) \approx \left(\left.\frac{\delta \pi}{\delta \zeta}\right|_{\zeta = \mu}\right)^2 \times \text{Var}\left(\zeta\right) = \left[\phi\left(\frac{\mu - q_{\alpha} \sigma}{\sqrt{\omega^2 + (1 - \rho)\sigma^2}}\right)\right]^2 \times \frac{\tau^2 + \rho \sigma^2}{\omega^2 + (1 - \rho)\sigma^2}. \] Thus, \[ \begin{aligned} \text{Var}(N_A) \approx m \times \psi \left(1 - \psi\right) + m (m - 1) \times \left[\phi\left(\frac{\mu - q_{\alpha} \sigma}{\sqrt{\omega^2 + (1 - \rho)\sigma^2}}\right)\right]^2 \times \frac{\tau^2 + \rho \sigma^2}{\omega^2 + (1 - \rho)\sigma^2}. \end{aligned} \] If the amount of common variation is small, so $\tau^2$ is near zero and $\rho$ is near zero, then the contribution of the second term will be small, and $N_A$ will act more or less like a binomial random variable with size $m$ and probability $\psi$. On the other hand, if the amount of independent variation in the effect sizes is small, so $\omega^2$ is near zero and $\rho$ is near 1, then the term on the right will approach $m(m - 1)\psi(1 - \psi)$ and $\text{Var}\left(N_A\right)$ will approach $m^2 \psi(1 - \psi)$, or the variance of $m$ times a single Bernoulli variate. So you could say that $N_A$ has anywhere between $1$ and $m$ variate’s worth of information in it, depending on the degree of correlation between the effect size estimates.

Interactive distribution

Here is an interactive graph of the probability mass function of $N_A$, with probability points calculated using Gaussian quadrature. Below the graph, I also report $\psi$, the exact mean and variance of $N_A$, and the first-order approximation to the variance (denoted $V_{approx}). When $\tau > 0$ and $\rho > 0$, the approximate variance is not all that accurate because the first-order approximation to $\text{Var}(\pi)$ isn’t that good.

Code

math = require("mathjs")
norm = import('https://unpkg.com/norm-dist@3.1.0/index.js?module')

quad_points = JSON.parse(all_quad_points).at(qp - 1)

sigma = 2 / math.sqrt(ESS)

zeta_sd = math.sqrt(tau**2 + rho * sigma**2)
ID_sd = math.sqrt(omega**2 + (1 - rho) * sigma**2)

crit = norm.icdf(1 - alpha)

binomial_coefs = Array(m+1).fill(null).map((x,index) => {
  return math.combinations(m, index);
})

probs = quad_points.map(zeta => {
  let Z = (zeta[0] * zeta_sd + mu - crit * sigma) / ID_sd;
  return [norm.cdf(Z), zeta[1]];
})

p_binom_norm = binomial_coefs.map((coef, a) => {
  let p = probs.map((x) => {
    return (x[0]**a) * ((1 - x[0])**(m - a)) * x[1];
  });
  return coef * math.sum(p);
})

Code

psi = norm.cdf((mu - crit * sigma) / math.sqrt(tau**2 + omega**2 + sigma**2))
psi_print = psi.toFixed(3)

E_NA = m * psi 
E_NA_print = E_NA.toFixed(3)

dpi_dzeta = norm.pdf((mu - crit * sigma) / ID_sd)
V_pi_approx = (dpi_dzeta * zeta_sd / ID_sd)**2
V_approx = m * psi * (1 - psi) + m * (m - 1) * V_pi_approx
V_approx_print = V_approx.toFixed(3)

V_NA = {
  let V_NA = 0 - E_NA**2;
  for (let i = 0; i <= m; i++) {
    V_NA += i**2 * p_binom_norm[i];
  }
  return V_NA;
}
V_NA_print = V_NA.toFixed(3)

Distribution of $N_A$

Code

Plot.plot({
  x: {
    label: "Number of significant effect sizes"
  },
  y: {
    domain: [0, 1],
    label: "Probability"
  },
  marks: [
    Plot.ruleY(0),
    Plot.barY(p_binom_norm, {
      fill: "steelblue"
    }),
  ]
})

Moments of $N_A$

Code

viewof m = Inputs.range(
  [1, 30], 
  {value: 6, step: 1, label: "Number of effect sizes (m):"}
)

viewof ESS = Inputs.range(
  [4, 300], 
  {value: 80, step: 1, label: "Effective sample size:"}
)

viewof mu = Inputs.range(
  [-2, 2], 
  {value: 0.3, step: 0.01, label: "Average effect size (mu):"}
)

viewof tau = Inputs.range(
  [0, 1], 
  {value: 0.1, step: 0.01, label: "Between-study SD (tau):"}
)

viewof omega = Inputs.range(
  [0, 1], 
  {value: 0.1, step: 0.01, label: "Within-study SD (omega):"}
)

viewof rho = Inputs.range(
  [0, 1], 
  {value: 0.6, step: 0.01, label: "Sampling error correlation (rho):"}
)

viewof alpha = Inputs.range(
  [0.005, 0.995], 
  {value: 0.025, step: .005, label: "One-sided significance threshold (alpha):"}
)

viewof qp = Inputs.range(
  [1, 30], 
  {value: 21, step: 1, label: "Number of quadrature points:"}
)

Back to top

Distribution of the number of significant effect sizes

Compound symmetry to the rescue

Just the moments, please

Interactive distribution

Distribution of \(N_A\)

Moments of \(N_A\)